Quadmesh

Michael D. Sumner

2018-08-08

Source: vignettes/quadmesh.Rmd

quadmesh

Build a “quadmesh”" in R.

## Loading required package: sp
## Warning in rgl.init(initValue, onlyNULL): RGL: unable to open X11 display
## Warning: 'rgl_init' failed, running with rgl.useNULL = TRUE

A “quadmesh” is a dense mesh describing a topologically continuous surface of 4-corner primitives. I.e. a grid, without the “regular”. This is useful particularly when combined with map projections and texture mapping.

We are not limited to a regular grid, trivially let’s distort the mesh by a weird scaling factor.

The topology of the grid is still sound, but we are no longer bound to the regular constraint.

## null 
##    3
shade3d(qm, col = grey(scl(qm$vb[3,qm$ib])))

rglwidget()

Why quadmesh?

Why meshes at all?

The simplest kind of mesh is a basic raster. Consider a matrix of 12 values.

##      [,1] [,2] [,3] [,4]
## [1,]    1    4    7   10
## [2,]    2    5    8   11
## [3,]    3    6    9   12

On its own this matrix has absolutely nothing to do with spatial data, it is literally a collection of 12 numeric values in a given order, and by the magic of programming we’ve nominated a shape of 3x4. We can’t help but think about this shape spatially however, but there’s a problem. Does each element occupy space or should we consider them to be infinitesimal locations?

R provides either interpretation (to simplify this story we nominate locations for the rows and columns explicitly).

When considered as an image, each matrix element occupies a certain space in width and height, but when considered as a point set the numbers simply float at the given locations. Which is correct? (Spoiler: Both are correct, it simply depends what we are doing.)

x <- seq(1, nrow(m)) - 0.5
y <- seq(1, ncol(m)) - 0.5
image(x, y, m)

text(expand.grid(x, y), lab = m[])

The raster package defaults to the image interpretation and helpfully assumes the values are nominally at the centre points as shown above. We have to nominate the extent or we end up in 0,1 range, we also have to invert the order of the values because raster counts from the top of the page and R’s matrix uses column-major order.

library(raster)
(r <- raster(t(m[, ncol(m):1]), xmn = 0, xmx =ncol(m), ymn = 0, ymx = nrow(m)))
## class       : RasterLayer 
## dimensions  : 4, 3, 12  (nrow, ncol, ncell)
## resolution  : 1.333333, 0.75  (x, y)
## extent      : 0, 4, 0, 3  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## data source : in memory
## names       : layer 
## values      : 1, 12  (min, max)

R’s image and rasters in general are so efficient because they only store this minimal amount of information: the actual data values, and the extent and dimensions of the space they occur in. If we had to store the centre coordinate of every cell, or worse the corner coordinates then the data storage goes up dramatically. Every software that deals well with these kinds of data has to treat these coordinates as implicit. We can easily expand the centre coordinates.

xyz <- as.data.frame(r, xy = TRUE)
head(xyz)
##           x     y layer
## 1 0.6666667 2.625    10
## 2 2.0000000 2.625    11
## 3 3.3333333 2.625    12
## 4 0.6666667 1.875     7
## 5 2.0000000 1.875     8
## 6 3.3333333 1.875     9
tail(xyz)
##            x     y layer
## 7  0.6666667 1.125     4
## 8  2.0000000 1.125     5
## 9  3.3333333 1.125     6
## 10 0.6666667 0.375     1
## 11 2.0000000 0.375     2
## 12 3.3333333 0.375     3

but to expand the corners we have to jump through some hoops and even then we get every instance of the corners, not only for each cell but to explicitly close the cell as a polygon.

## class       : SpatialPointsDataFrame 
## features    : 60 
## extent      : 0, 4, 0, 3  (xmin, xmax, ymin, ymax)
## coord. ref. : NA 
## variables   : 4
## names       : layer, Lines.NR, Lines.ID, Line.NR 
## min values  :     1,        1,        1,       1 
## max values  :    12,       12,        9,       1

There are only 20 unique coordinates at the corners, which is where quadmesh comes in.

## List of 6
##  $ vb           : num [1:4, 1:20] 1.49e-08 3.00 1.10e+01 1.00 1.33 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:4] "x" "y" "z" "1"
##   .. ..$ : NULL
##  $ ib           : int [1:4, 1:12] 1 2 6 5 2 3 7 6 3 4 ...
##  $ primitivetype: chr "quad"
##  $ material     : list()
##  $ normals      : NULL
##  $ texcoords    : NULL
##  - attr(*, "class")= chr [1:2] "mesh3d" "shape3d"

This is a mysterious seeming data structure, it is the mesh3d type of the ‘rgl’ package, rarely seen in the wild.

The structure is vb, the coordinates of the mesh - these are the actual corner coordinates from the input raster.

image(r)
op <- par(xpd = NA)
text(t(qm$vb), lab = 1:20)

Notice how these are unique coordinates, there’s no simple relationship between the cell and its value and its four corners. This is because they are shared between neighbouring cells. The relationship is stored in the ib array, this has four rows one for each corner of each cell. There are 12 cells and each has four coordinates from the shared vertex pool. The cells are defined in the order they occur in raster.

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,]    1    2    3    5    6    7    9   10   11    13    14    15
## [2,]    2    3    4    6    7    8   10   11   12    14    15    16
## [3,]    6    7    8   10   11   12   14   15   16    18    19    20
## [4,]    5    6    7    9   10   11   13   14   15    17    18    19

It works directly with rgl functions.

The primary means to create this format from a raster is for 3D plotting, but because we have access to the coordinate directly it provides other uses. We can transform the coordinates (i.e. a map projection) or manipulate them and augment the Z value (for example) in flexible ways.

(The usual way of driving rgl grid surfaces is rgl.surface but this is limited to the x, y, z list interface just as image() is. )